Base Question Implementation



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint



In [2]:

    
from initial_velocities import velocities_m, velocities_S
from DE_solver import derivs, equationsolver

Finally, I am ready to answer the base question that I have laid out, which is how a system of many stars would react when the disrupting galaxy came into a close orbit to the main galaxy. Toomre and Toomre's paper lays out the initial orbit of 120 stars around the main galaxy, 12 particles with a circular orbit of radius 20 percent $R_{min}(25\hspace{1 mm} kpc)$, 18 particles at 30 percent, 24 at 40 percent, 30 at 50 percent, and 36 and 60 percent.

Defining emtpy initial condition array:



In [3]:

    
ic_base = np.zeros(484)

Setting values for S, M, and t:



In [4]:

    
max_time_base = 1.5
time_step_base = 120
M_base = 1e11
S_base = 1e11



In [5]:

    
S_y_base = 70
S_x_base = -.01*S_y_base**2+25
vxS_base = velocities_S(M_base,S_base,S_x_base,S_y_base)[0]
vyS_base = velocities_S(M_base,S_base,S_x_base,S_y_base)[1]

Setting initial condition array values pertaining to S:



In [6]:

    
ic_base[0] = S_x_base
ic_base[1] = S_y_base
ic_base[2] = vxS_base
ic_base[3] = vyS_base

Creating positions for the stars of M:



In [7]:

    
particles = np.array([12,18,24,30,36])
percent = np.array([.20,.30,.40,.50,.60])



In [8]:

    
a=np.array([particles,percent])



In [9]:

    
x = []
y = []



In [10]:

    
for i in range(0,5):
    theta = np.arange(0,2*np.pi,(2*np.pi)/a[0,i])
    r = 25*a[1,i]
    for t in theta:
        x.append(r*np.cos(t))
        y.append(r*np.sin(t))



In [11]:

    
x_y = np.array([x,y])

The following plot shows all the initial positions of the stars:



In [12]:

    
plt.figure(figsize=(5,5))
for n in range(0,120):
    plt.scatter(x_y[0,n],x_y[1,n])

I also wrote these positions to disk so I could use them for other cases:



In [13]:

    
np.savez('star_positions.npz',x_y)

Putting these values into my initial condition array, as well calling the initial velocity function on each position:



In [14]:

    
for i in range(0,120):
    ic_base[(i+1)*4] = x_y[0][i]
    ic_base[((i+1)*4)+1] = x_y[1][i]



In [15]:

    
for n in range(1,int(len(ic_base)/4)):
    ic_base[n*4+2] = velocities_m(M_base,ic_base[n*4],ic_base[n*4+1])[0]
    ic_base[n*4+3] = velocities_m(M_base,ic_base[n*4],ic_base[n*4+1])[1]

Calling my differential equation solver, and saving the data to disk:

Times vary, either just under or just over a minute



In [39]:

    
sol_base = equationsolver(ic_base,max_time_base,time_step_base,M_base,S_base)



In [40]:

    
np.savez('base_question_data.npz',sol_base,ic_base)



In [ ]: